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I. INTRODUCTION 


The trend in packaging of electronic components is to 
place more power in an ever decreasing space. Computers that 
occupied an entire room 20 years ago, now occupy a small 
corner of that room. With the size of chips decreasing, the 
volumetric power generation (q’”) is of the order of MW/m’ in 
magnitude. 

This can lead to high temperature operation of the chips 
with a resulting degradation in performance and or failure. 
The military has set the maximum temperature of a 
semiconductor junction to between 100 and 110 °C [{Ref. 1]. 

This investigation deals with the steady state temperature 
distribution and resulting thermally induced stresses in a 
tri-material electronic package as shown in Figure 1.1, 
dimensions not shown to scale. The package consists of a 
semi-conductor (chip) attached to a substrate material by a 


solder joint. The chip is a source of thermal energy. The 


package is exposed to a convective environment on all sides. 


¢ 


Ne Ly oes oT 


| SUBSTRATE 





Figure 1.1 Tri-material electronic package 


With the increase in volumetric power generation the 
overall temperature of the package is increased. The stresses 
in the package arise from two sources, (a) thermal strains 
(€,) caused by temperature gradients within a material and (b) 
differences in the materials mechanical properties, 
coefficient of thermal expansion (a) and stiffness (E). At 
material interfaces, displacement continuity requires 
expansion or contraction accommodations between materials. 
There are two such interfaces where this will occur, at the 


chip/solder interface, and at the solder/substrate interface. 


A. BACKGROUND 

Thermal stresses of bi-material assemblies were 
investigated by Timoshenko in a 1925 paper [Ref. 2], where a 
bi-metallic thermostat assembly was subjected to a uniform 
temperature field. Timoshenko obtained analytical expressions 
for the bending stresses through the cross section of the 
bimetallic strip. He also noted that the ‘distribution of 
shearing stresses along the bearing surface cannot be 
determined in an elementary way, that they are of a local type 
concentrated near the ends of the strip’, and that the 
shearing stresses can be of the same magnitude of the bending 
stresses. He mentions, in passing, the existence of ‘local’ 
normal stresses between the assembly interfaces but says 
nothing of their magnitude. Goland and Reissner [Ref. 3] 


determined the shear and peeling stresses at the interface 


of a cemented lap joint, where the peeling stress tends to 
pull the materials apart. 

In a comprehensive review of the subject of tri-material 
assemblies, Suhir [Ref. 4] took the analysis one step further 
by investigating the stresses in a tri-material system due to 
a thermal environment. His model is based on an assembly 
fabricated at elevated temperatures and subsequently cooled. 
Specific restrictions were placed on the tri-material assembly 
in a uniform temperature field. He obtained analytical 
solutions of the problem. For the present study the general 
tri-material problem is solved numerically. In particular the 


tri-material system is that of an electronic package. 


B. PROBLEM DEFINITION 

The problem was partitioned into two parts, a thermal 
analysis, whereby the heat conduction equation was solved, and 
a thermal stress analysis whereby the stresses produced by the 
thermal field were obtained. Both problems were solved 
numerically by the Finite Element Method (FEM). For the 
numerical analysis, symmetry is invoked along the centerline 
of the system as shown in Figure 1.1. The steady state 
thermal study focuses on the effects material properties, 
geometric configuration, convective cooling, and encapsulation 
have on temperature distribution in an electronic package. 
Once the temperature field is known the data is input into a 


program that assigns temperatures to another FEM mesh which is 


used in the stress analysis. The stresses examined include 
bending stress at the centerline of the package (9o,), the 
bending (om) normal (o,,), and shear (t,), stresses along the 
chip/solder interface, and the bending (o,,), normal (o,,), and 


shear (t,) stresses at the solder/substrate interface. 


II. FEM DEVELOPMENT OF THERMAL PROBLEM 


A. PROBLEM STATEMENT AND ASSUMPTIONS 

As noted in the introduction, the purposes of this 
investigation are to (a) determine the variation of 
temperature in the electronic package and (b) determine the 
stresses associated with these temperatures. The results here 
are the steady state results. 

The steady state thermal study focuses on effects that 
certain key parameters have on system behavior. They include, 
the effects total power (P) in the chip, material 
combinations, convective cooling, geometric configuration and 
encapsulation have on the temperature distribution in the 
electronic system. 

Finite element codes were used to solve for the 
temperature distribution and stresses of the electronic 
package. In order to limit the mathematical complexity of the 
problem the following approximations and specific assumptions 
for FEM program HEATSTEADY in Appendix A. are made as follows: 

¢ The thermal contact resistances due to ‘imperfect’ contact 
at the solder interfaces are negligible. 
¢ The thermal conductivity and coefficient of thermal 


expansion (k and a) are considered constant over the range 
of temperatures encountered in the study. 


- Because the temperatures are not very large radiation heat 
transfer between the package and the surrounding enclosure 
is negligible. 

- The convection heat transfer coefficient (h) is constant 
over the entire edge upon which it acts. 

All of these assumptions are reasonable within the scope 
Cfecthasmistudy: Once the temperature distribution of the 
package is obtained the position and temperature of all the 
Global nodal points are entered into program Filter in 
Appendix B. to obtain temperatures that correspond to FEM 
program Weld in Appendix C. to determine the stresses 


developed in the package. The development of program WELD is 


discussed in Chapter III. 


B. FEM FORMULATION OF THERMAL PROBLEM 

The Galerkin FEM 1S an approximation method which 
transforms a linear partial differential equation into a 
system of linear algebraic equations. Using the two- 


dimensional heat equation [Ref. 5]: 


V- (kV 7T) +q"=0 (2.1) 


and boundary conditions 


Cauchy: ~K OT =h(T-f,) (2.2) 
dn 
Nuemann: - 0 (2.3) 
_ on 


where the V operator denotes wae By , and thus 
Vi(kKVT) = kVT= OT , ST (2.4) 
Ox Oy 


where k is the thermal conductivity, T is the temperature, and 
q’’ is the volumetric heat generation. The temperature 
distribution can be determined at system nodal points of the 
electronic package. For this method, a linear triangular 
shape function (N,;) which possesses the Kronecker Delta 


property (equation 2.5): 


N,;(nodej) =6,, 1£i1=]j7 then N,;=1 


Ji Fi 2 
ame: #9 then N, = 0 Med 


is used. The linear shape functions maintain function 
continuity throughout the domain. 

Heat flux continuity at material interfaces is built into 
the FEM formulation. The thermal FEM grid and boundary 
conditions are given in Figure 2.1 and Figure 2.2 


respectively. 
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B DRY BOUNDARY CONDITION 





BOUNDARY CONDITIONS FOR 
THERMAL PROBLEM 


Figure 2.2 FEM thermal problem boundary conditions 


An approximate solution, t, for temperature, T(x,y), as 


formed as follows: 


B= 6 Nios (2.6) 


R 


where T is the exact solution of the heat equation in 
continuous space, t is the approximate solution in discrete 
space, {N}' is the transpose of a column vector of the 
triangular linear shape functions, and {T} is the vector of 
nodal temperatures. 

After the approximation is formulated, the next step is to 


form the Residual, R, as follows: 


R=2(t) - q(s) (2.7) 


where q’’ is the heat generation term and Y& denotes the 
differential operator which in the case of the heat equation 


is defined by: 


9(t) = V(kVE) (2.8) 


With this substitution, the residual becomes: 


R= (V-(kKV (Ni {7t) )] - q(s) (2.9) 


From the residual, the Galerkin Equations are formed: 


1 {w}(R) ds = {0} (2.10) 


10 


where D denotes surface integration and where {0} is the null 
vector. Further substitution for R into the Galerkin vector 


equation results in: 
[im (V-(kV ({w}T{ TH) ) ds - [ ima (s) ds = {0} (2.11) 
D 


To solve the Galerkin Equation, Green’s Theorem is evoked 


which yields: 
¢ {w) CATE els ee { [(Vin} (kV (imi {T}) )] ds + [ (wtads= {0} 
B D D 


where B subscript on the integral denotes evaluation of these 
integrals around the boundary of domain D of the electronic 
package. When integrated the boundary integral and the heat 
generation integral form the “force" vector {F}. The middle 
integral forms the [A] matrix. Thus the heat conduction 


differential equation becomes, 


[A] ir} = {PF} (2.13) 


where {T} contains the vector of nodal temperatures. 


C. THERMAL PROGRAM DESCRIPTION AND VALIDATION 

The main thermal program, HEAT.FOR, a VAX Fortran 77 Code 
for the FEM thermal analysis is contained in Appendix A. The 
program begins by reading in a data file DATA3.DAT. This 
input data is then processed through several subroutines that 


perform the FEM analysis. 


ee 


1. Subroutine Input Description 
The following information is supplied to the Input 
subroutine from a data file. 


Input file form: NCOL1,NCOL2,NROW, 
XPOS, YPOS, DIVX, DIVY, 
THERCON,GEN,QTPR, 
DIVSUB, DIVXAIR, DIVYAIR, 
NBCLOW, NBCUP, NBCVERT, NBCSUB, TAMB, 


HLOW, HUP, HVERT, HSUB 
Parameter Description 
NCOL1 The number of columns of substrate in the FEM grid. 
NCOL2 The number of columns of chip in the FEM grid. 
NROW The number of rows in the FEM grid. 
XPOS Starting position for new grid density along xX 
axis. 
YPOS Starting position for new grid density along Y 
axis. 
DIVX Spacing of grid density along X axis. 
DIVY Spacing of grid density along Y axis. 
THERCON Thermal conductivity (k). 
GEN 0 no heat generation term 


Hol 


1 heat generation term 
OTPR Amount of volumetric heat generation (qd). 
DIVSUB Number of rows in FEM grid composed of substrate. 
DIVXAIR Number of columns in FEM grid composed of air. 
DIVYAIR Number of rows in FEM grid composed of air. 
NBCLOW O = insulated boundary at the bottom of substrate. 
1 = convective boundary at the bottom of substrate. 
NBCUP 0 insulated boundary at the top sof chap: 


l = convective boundary at the top of chip. 
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insulated boundary on the side of the chip. 


1 convective boundary on the side of the chip. 
NBCSUB 0 = insulated boundary on top of the substrate. 

1 = convective boundary on top of the substrate. 
HLOW Convection coefficient corresponding to NBCLOW. 
HUP Convection coefficient corresponding to NBCUP. 
HVERT Convection coefficient corresponding to NBCVERT. 
HSUB Convection coefficient corresponding to NBCSUB. 


2. Subroutine Grid Description 

Spleen UNE GRID takes the input information and 
constructs a grid of horizontal and vertical lines on a 
geometry as shown in Figure 2.1. This grid will be used to 
form the triangular elements used in the thermal FEM analysis. 
This subroutine allows the user to generate a variety of 
meshes and/or refinements and permits validation of grid 
independence very quickly. 

To optimize computer time the subroutine is able to 
generate different grid meshes with a variety of spacing 
options as shown in Figure 2.3. Fine meshes are used at 
convective boundaries and material interfaces. A coarse mesh 
is used within a material away from boundaries. The 
Subroutine allows the user to independently vary the mesh 
density both horizontally and vertically for up to nine 
different densities in either direction. The number of mesh 
densities can be increased by simply adding a number of "IF 


THEN" statements at the top of the program (line 34). 
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fine coarse coarse 


fine 


Figure 2.3 Fine and coarse vertical mesh and fine and coarse 
horizontal mesh 


Local nodal points are assigned to each element, one 


at each corner, counterclockwise as shown in Figure 2.4. 





Figure 2.4 Triangular element numbering 

Next the subroutine generates a correspondence table 
which identifies the correspondence between local nodal points 
and global nodal points. Next, each triangular element is 
identified with a specific material and assigned a thermal 
conductivity (k) and whether or not it will have a volumetric 
heat generation term (q’”’). Lastly boundaries are assigned 


according to Figure 2.2. 
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3. Subroutine Heatmat, Heatmod and Lsarg Description 
Subroutine HEATMAT constructs the global matrix [A] 
and the heat generation vector {F} according to the 
correspondence table made in subroutine GRID. 

Subroutine HEATMOD is not used in this program but can be 
used to modify the global matrix [{A] for a Dirichlet boundary. 
HEATMOD replaces the Galerkin equation for the specific nodal 
points with the specified temperatures. LSARG is an equation 
solver provided by the IMSL Library used to solve equation 
2el3. ° 

4. Subroutine Output Description 

Subroutine OUTPUT creates data files and performs an 
energy balance between the thermal energy which is convected 
at the boundaries and the thermal energy generated by the chip 


according to: 


4 
Pee a, = O° Aca, (2.14) 
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where h, and 1, correspond to the convection coefficient and 
length of each convective boundary, T, corresponds to the 
average temperature across the boundary and A,,;, 1s the area 
of the chip. This heat balance is conducted for each case 


study. 
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D. VALIDATION I. DIRICHLET AND CAUCHY BOUNDARY 

A two dimensional body subjected to both Dirichlet and 
Cauchy boundary conditions was taken from [Ref. 6] and run on 
program HEAT to validate the programs ability to handle both 
boundary conditions. The results of program HEAT were 


identical to Reference 6. 


E. VALIDATION II. GRID INDEPENDENCE 

In order to ensure that the numerical errors are kept to 
a minimum the size of the mesh was decreased with incremental 
steps until there was a less than 1% change in the solution 
from one grid to the next finest grid (grid independence) . 

The final mesh chosen has 125 degrees of freedom (DOF) 
with 202 elements. The mesh which provided less than 1% 
difference had 161 DOF with 268 elements. The total power of 
the chip used is P = 40 W. The total power convected for the 
DOF = 125 grid was P,,,, = 39.95 W. The total power convected 


for the DOF = 161 grid was P 39.94 W. The percent 


conv 


difference in total power P from the coarse mesh to the fine 


mesh is 0.03%. 


F. VALIDATION III. ENERGY BALANCE 

As mentioned previously, there is an energy balance in 
subroutine OUTPUT which is given by equation 2.14. This 
checks program HEAT’s ability to correctly handle a generation 


term. The total power for the majority of case studies was 
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equal to P = 1W. The average convective power output was 


= .99785 W for a percent difference of 0.21%. 


qeonv 
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III. FEM STRESS DEVELOPMENT 


The finite element stress formulation utilizes a recently 
developed element which provides for axial and lateral 
displacement continuity. The results of the stress code are in 
good agreement with existing solutions [{Ref. 2,4] to problems 


with uniform temperature fields. 


A. FORMULATION OF THE MODEL 

A brief description of the FEM formulation for stresses 
follow. In Figure 3.1, each element has six degrees of 
freedom; axial displacements at the four corner points, and 
lateral displacements at the two ends. An advantage of the 
element is that axial and lateral displacement continuity 


results. 


v, — > ms 
Ni wif 


Figure 3.1 A typical element with 6 degrees of freedom 
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The axial displacement field u(x,y) is assumed to be linear in 


both the axial and transverse directions. That is, 


2 
u(x,y) = D.N, (x)[H, (y) up +H, (y) uz] (3.1) 


1=1 


Superscripts b and t on the nodal displacements refer to 
bottom and top displacements respectively. The linear shape 


functions N, and H, are given by 


H, (y) -a¥ H, he N, (x) = =—* N, == 


The lateral displacement field is given by 


N, (X) Vv; (3.2) 


Me 


v(x) 


The strain-displacement relations are, 














H. H. 
€, = se = te [ug - uy] + ae Ex = uy" (3.3) 
and 
> N. (x) b N. (x) b Vine ie 
Mosy @ a Ee ee) = a —— (3.4) 
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Defining the axial displacement vector as 


{8 }7 = ( ee us’ Ue ce ) (3.5) 


and the stiffness matrix as 


ix] = ff “eh etal dy dx (3.6) 


where the B vector is, 


ON. ON. ON. ON. 
7° i 1 2 2 3.7 
{B}? = ( erp oh H,) (3.7) 


and the force vector due to temperature as, 


_ Ehpira 
(F,} = aa i {B} Fa aTdy dx (3.8) 


gives the bending matrix equations as, 


[K,] {6,} = Fj (3.9) 


Equation (3.9) defines the bending behavior. Behavior due to 
shear is obtained as follows. Define the row vector of 


displacement degrees of freedom as 


SJ? =(up uy vv, Sa eee (3.10) 


The shear stiffness matrix 1s given by 


20 


(K.] = [of (ede ay ax (any 
PA o Jo 


where 


OH, dH, ON, dH, dH, 


i a o_—e — es as 
eee G25) 82 Sy 


ON2) (3.12) 
Ox 


which gives the equations for shear behavior 


[K,] {6,} = {o} (2513) 


The matrix equations for bending and shear behavior are 


combined to give the stiffness equations for the system. 


B. PRELIMINARY RESULTS 

Assessment of the present model was obtained from 
comparisons with previous analytical models (Ref. 2,4]. Exact 
agreement for bending stresses was obtained when the present 
model was applied to a sample bimetallic case solved in 
[Ref.2]. For a tri-metallic case the present model produced 
shear stress results in reasonably good agreement with results 


obtained from [Ref. 4]. 
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IV. STEADY STATE TEMPERATURE ANALYSIS 


A. OVERVIEW OF STEADY STATE TEMPERATURE ANALYSIS 

The analysis begins with an FEM approximation of the 
temperature distribution of the tri-material electronic 
package. First the chip volumetric power generation is 
examined to determine the effect an increase in q’”’ has on the 
chip temperature. Keeping the power of the chip and the level 
of convective cooling constant, several packaging materials 
are combined to determine how the temperature distribution 
varies with the variation in the material thermal 
conductivities. This is followed by a convection cooling 
analysis, in which temperatures and temperature gradients are 
examined for a common electronic package while the convection 
coefficient h varies. The package geometry is then studied, 
holding the overall package height equal and varying the chip 
size as a percentage of overall package height. This is done 
for the case of constant heat flux, and for the case of 
constant volumetric heat generation. Lastly the effects of 
placing a protective coating "encapsulation" on the chip is 
studied to determine the variation in temperature distribution 


from a chip which has not been encapsulated. 
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B. EFFECT OF VOLUMETRIC POWER GENERATION ON TEMPERATURE OF 

THE CHIP 

The heat fluxes of electronic packages is expected to 
reach nearly 10° w/m? by the year 2000. This is two magnitudes 
in power over the chip studied here. The total power output 
(P) of the chip studied is one watt. This translates to a 
volumetric power generation (q’”) of 20 MW/m’ for a device 5mm 
x 5mm x 2mm in dimension or a heat flux (q’) of 4x10* W/m*. 
The dimensions chosen are of a typical electronic device size. 

The temperature results throughout the study can be 
modified for any power output due to the linearity of the 
generation term in the heat equation. In Figure 4.1 the 
temperature of the chip is plotted against chip volumetric 
power output (q’’’). All other quantities must be constant 
for this to remain true, including boundary conditions and 
other materials within the package. For the particular 


geometry studied, temperature (T) as a function of volumetric 


power output q’” W/m’ is given by equation 4.1 


T=m*q"@+ T. (4.1) 


where the change in temperature per unit change in volumetric 
power output is b = 4.65 °C /MW/m?* and T, 1s the ambient 


temperature of the convective fluid. 
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Figure 4.1 Chip temperature as a function of volumetric power 


C. EFFECT OF VARIOUS MATERIAL COMBINATIONS 


The materials used in this investigation along with their 
key thermal properties are listed in Table 4.1. The set of 
materials that comprise the electronic tri-material package 
focused on for the majority of this study are the silicon 
chip, a Pb-Sn solder, along with several substrate materials 


including, Epoxy Fiberglass, Polyimide Fiberglass, Alumina, 
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and Aluminum Nitride. These are the most common materials 
used in electronic packaging. The other chip and solder 
materials considered in this study account for _ large 
variations in thermal conductivity and/or thermal expansion 
coefficient (a). 


TABLE 4.1 ELECTRONIC PACKAGE MATERIALS AND PROPERTIES 


Use material k (W/mK) opel OEY aly 

Chip Silweon 150 Ze 
Gallium Arsenide 58 57 

Substrate Epoxy Fiberglass 0.16 140 - 180 
Polyimide Fiberglass G35 P20 Go 
Alumina 18 60 
Aluminum Nitride 230 33 

Solder Pb(5%)/Sn 63 290 
Au(20%)/Sn | 57 159 

Air Gap 0.026 - 


Throughout the study within this section the following 
parameters are held constant; volumetric power (q’’) of the 
chip, ambient temperature T, of the convective fluid, and 
coefficient of convection cooling (h). Figure 4.2 shows the 
non-dimensional temperature profiles of the package at a non- 
dimensional distance, ¢, = 0.564. The plotted values of T, 
were obtained at the FEM grid nodes. The temperatures which 
lie along ¢, = .564 were chosen because they pass through the 
solder joint of the package. Thermal stresses will result not 
only from temperature gradients within a material, but also 
from the direct contact of the three different materials 
within the package which have large differences in their 


coefficients of thermal expansion (a). In Figure 4.2 the 
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abscissa is a normalized temperature T, given by equation 4.2, 


—— (4.2) 


~ 


where T, is the FEM grid temperature. The ordinate is he 


non-dimensional y location ¢, that 1s given by equation 4.3. 


oy Se (4.3) 


where H is the total height of the package and y is the FEM 


grid global nodal location. 


1.Q 
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Figure 4.2 Non-dimensional temperature for silicon chip, Pb-Sn 
solder and various substrates. 
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All sixteen combinations of materials were analyzed fora 
temperature profile. There was no significant difference in 
temperature profile with regard to which chip or solder 
material was used in the package. For the volumetric power 
level used in this analysis the difference in thermal 
conductivity of the Silicon chip (k=150 W/m°C) or Ga/As chip 
(k=58 W/m°C) does not influence the temperature profile. A 
study was conducted in which the volumetric power was raised 
by four times the initial amount. This failed to produce any 
Significant difference in temperature profile or magnitude. 

The temperature profiles show that the substrate is the 
governing material which determines the temperature profile 
a given Biecereniic package takes on. The following symbology 
will be used in further analysis. A subscript of 1 refers to 
the chip, subscript 2 refers to the solder, and subscript 3 
refers to the substrate. Low thermal conductivity substrates 
used in this analysis have ratios of k,/k, < 0.0025. High 
thermal conductivity substrates used in this analysis have a 
ratio range of 0.10 < k,/k, < 1.6. 

Meciatacteristic Of a high k,/k, package is a uniform 
temperature field throughout all materials in the package. 
All materials within the package can conduct heat very well. 
With the low k,/k, substrates a temperature gradient is 
developed within the substrate in both the X and Y directions. 
As a result of not being able to conduct the heat effectively 


through the low conductive substrate, the chip temperature is 
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increased by approximately 40% over a package with a highly 
conductive substrate. 

The solders have thermal conductivities which are close to 
the chip thermal conductivity (k,/k, = 0.446) and do not effect 
the temperature profile. However, the coefficient of thermal 
expansion a between the two solders chosen vary by 55% thereby 
influencing the thermal strain profile significantly. The 
solder joint takes on the chip temperature in all cases 
because of its ability to conduct the heat generated by the 
chip. 

Thermal strain profiles are presented in Figures 4.3 and 


4.4 for the two different solders. In these figures, the 
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Figure 4.3 Thermal strain profile for silicon chip, Pb-8Sn 
solder and various substrates. 
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abscissa is the difference of the nodal temperature and 
initial temperature multiplied by the coefficient of thermal 


expansion a given by 
ee, 5 peo (4.4) 


and the ordinate is given by equation 4.3. 


¢ 20 MW/m? 
= 250 W/m? °C 
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Figure 4.4 Thermal strain profile for silicon chip, Au-Sn 
solder and various substrates. 
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1. Low Thermal Conductivity Substrate Packages 

Further analysis is conducted by dividing the packages 
into two groups, low thermal conductivity substrates and high 
thermal conductivity substrates. The two substrates chosen, 
Epoxy Fiberglass (k,/k, = 0.00106) and Polyimide Fiberglass 
(K,/k, = 0.0023), are part of a group of substrates that have 
k, ranging from 0.12 W/m'C to 0.35 W/m C [Ref. 1]. The 
temperature achieved throughout the chip for each of the case 
studies, Figure 4.2, are approximately constant, that is, not 
dependent on position. The variation in temperature of the 
chip was limited to less than 1.2% over the entire cross 
section. The solder achieved the same temperature as the chip 
due to the solders high thermal conductivity. The variation 
between the chip temperature and the solder temperature is 
less then 0.1%. 

Within the substrate temperature gradients 9, and 9, 
are developed in both the X and Y directions respectively. 


Temperature gradients are given by: 


oye! es i aca! 
Ox 7 8. (xX, = so) (2a 
oF tg . “477 (4.6) 


oy (Yj eee 


where T; and T, are FEM Global nodal temperatures, x, and Ma 
are FEM Global nodal distances in the X and Y directions 


respectively. 
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These temperature gradients are only found in low 
thermal conductivity substrates and an explanation for this 


can be found by examining Fourier’s Law, equation 4.7, 


of! = mee = ke 


— (4.7) 


x 


where q’’ is the heat flux in the X direction. 

The solder acts as a sort of ‘heat funnel’ for the 
chip to transfer its thermal energy to the substrate since the 
air pocket under the chip acts as an insulator. By examining 
Fourier’s law there are two ways to increase the heat flux q”. 

First by increasing the thermal conductivity (k) of 
the material, and secondly by increasing the temperature 
difference (or temperature potential). The convection heat 
transfer from the chip surface to the surrounding medium is 
unable to transfer all the thermal energy in the chip 
generated by the power source q’’’ and therefore results in the 
conduction of heat flux through the solder. When the thermal 
conductivity of the substrate is low (as in the Fiberglass 
substrates), the temperature difference through the substrate 
must become proportionally larger to account for the heat flux 
out of the chip. For this reason, when the substrate has a 
high thermal conductivity, the temperature potential need only 
be very small to conduct the same amount of thermal energy. 

In the X direction the gradient is largest at the 


upper right hand side of the solder/substrate interface. As 


Syl 


Figure 4.5 shows, the gradient 9, diminishes by roughly 72% 
when only 25% into the substrate. The temperature gradient 
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Figure 4.5 6, in epoxy fiberglass substrate 

has decreased by 88% when half way through the substrate. 
This behavior can be accounted for by the fact the thermal 
resistance of the substrate in conduction is greater then the 
thermal resistance of the convection surface at the top of the 
Substrate. Equivalent thermal resistances in conduction and 


convection are given by: 


Bie 


= —— (4.8) 


where Ax is the distance between two points of interest 


R (4.9) 


9 ee! 
conv — RA 


and A is the area in which the heat flux is acting. 

Taking the area under the solder joint for conduction 
and the top of the substrate to the right of the solder joint 
for convection the following equivalent thermal resistances 


are obtained: 


0.002 ° 
R ee 14.2 °C/W. 4.10 
cence OmiG(G. 00088) (1) / ( 


di: 


R i 
sch BOO (0,001) (1) 


= 5 °C/W (4.11) 


which shows clearly that the thermal energy has approximately 
1/3 less resistance at the convective boundary versus 
conduction through the substrate. 

The gradient in the X direction on the left hand side 
of the solder joint is only 41% of that on the right, as shown 
in Figure 4.6. This again is a result of a higher thermal 
resistance in conduction within the substrate. With the 
insulated boundary condition imposed in the FEM formulation, 
the gradient in the X direction at X=0 and X=L is necessarily 


zero. 
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Figure 4.6 6, in epoxy fiberglass substrate 


The temperature gradient in the Y direction, O., is 
largest at 49.02 ‘C/mm at the solder/substrate interface as 
shown in Figure 4.7. Away from this interface the gradient 
decreases to approximately 20 ‘C/mm and remains constant 


throughout the entire substrate. 
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Figure 4.7 6, in epoxy fiberglass substrate 


The position where the largest temperature gradients 


in both the X and Y directions are found are shown in Figure 


4. 


8. 


These are the points where maximum heat fluxes are 


conducted through the solder/substrate interface. 
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Figure 4.8 Maximum temperature gradients in epoxy fiberglass 
substrate 
2. High Thermal Conductivity Substrate Packages 
The two substrates chosen in this range are Alumina k, 
= 18 W/m'C (k,/k, = 0.12) and Aluminum Nitride (ALN) k, = 230 
W/m'C (k,/k, = 1.53). Referring back to Pigure 4.2 themieis 
dimensional temperature distribution along ¢, = .564 for 
Alumina and AI1N graph is approximately a vertical straight 
line. Therefore there are no temperature gradients anywhere 


within the package. The reason for no temperature gradients 
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in the substrate is due to the substrate’s high thermal 
conductivity, enabling it to diffuse a large heat flux without 
developing a temperature gradient, as shown by Fourier’s Law 
of heat conduction, equation 4.7. Figure 4.3 shows the 
largest shear strains at the solder interface have been 
reduced by 40% over the low thermal conductivity substrates 
due to the increased ability of the high thermal conductivity 
substrates to conduct heat away from the chip thereby lowering 
the overall temperature of the package and decreasing the 
thermal expansion. 

There is however no appreciable temperature difference 
between the Alumina and AI1N substrate despite the thermal 
conductivity of the Aluminum Nitride being approximately 13 
times as large as the Alumina thermal conductivity. The 
limiting factor for reduction of temperature of the package 
has become the convective thermal resistance, which must be 
lowered in order to lower the package temperature. The 
conductive thermal resistances of the Alumina and AIN 


substrate and convective thermal resistance follow: 


0.002 A 
Alumina R ie eee a 0), 13° .C/W 
cond = s.18 (0.00088) (1) / 


0.002 P 
AIN R eee = ) 1 C/W 
Soho” Saepaig"( 0, 0008 Rome ) / 


Si 


R ee 
ete 2 0.0s60..00)) cle) 


= 5°C/W 
showing that the thermal resistance to heat transfer by 


convection is 38 and 500 times as large as resistance to heat 


transfer by conduction in the Alumina and AI1N respectively. 


D. EFFECT OF CONVECTION HEAT TRANSFER ON SYSTEM BEHAVIOR 

This study focuses on the effect of convection cooling on 
the temperature gradients developed within the low thermal 
conductivity substrates and the overall temperature decrease 
of the “chip. The electronic packages with high thermal 
conductivity substrates do not develop any temperature 
gradients and exhibit an overall decrease in package 
temperature. The following parameters are held constant; 
volumetric power output (q’” W/m’) of the chip, ambient 
temperature (T,) of the convective fluid, and the following 
materials, silicon chip, Lead/Tin solder and Epoxy Fiberglass 
substrate. 


The governing equation for convection cooling is given by: 


gat (Tr ar) (4.12) 


where q’” is the heat flux in W/m? and T is the surface 
temperature. The inverse relation of system temperature to 
convection coefficient is shown in Figure 4.9. This figure 
illustrates the need for ever increasing cooling to maintain 


electronic packages within a temperature tolerance range. An 
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increase in convective cooling from 100 W/m*K to 300 W/m°K 
results in a reduction of chip temperature of 90 ‘C. The same 
increase in convective cooling from 300 W/m*K to 600 W/mK 


produces a decrease in chip temperature of only 26 ‘C. 
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Figure 4.9 Convection coefficient vs. chip temperature 

The temperature gradients in both X and Y directions that 
are developed within the low thermal conductivity substrates 
are diminished by increasing the convective heat transfer as 


shown in Figures 4.10 and 4.11. 
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Figure 4.10 6, in epoxy fiberglass substrate vs. convection 
coefficient 
In the X direction the maximum temperature gradient is 
decreased by 60% by increasing the convective cooling by 250%. 
In the Y direction the maximum temperature gradient is 


decreased by 52% by increasing the convective cooling by 250%. 
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Figure 4.11 6, in epoxy fiberglass substrate vs. convection 
coefficient 


E. EFFECT OF GEOMETRY ON SYSTEM BEHAVIOR 

The dimensions of the chip and substrate are varied by 
assuming the chip is a percentage of a fixed overall height 
and the solder and substrate make up the remaining percentage. 
We define the non-dimensional term B as the ratio of chip 


height to overall package height (H) as in equation 4.13, 


41 


Ae iP 
= _ chip 4.13 
B H ( ) 


Substrates are approximately kept at a minimum of one half 
a millimeter (Ref. 1]. The following material combination was 
selected; a Silicon chip, a Lead/Tin solder and a Epoxy 
Fiberglass substrate. The following parameters were held 
constant; convective cooling, ambient temperature, the height 
of the air/solder interface and the overall height of the 
package. 

First, a study with a constant heat flux (q’’) was 
conducted by varying the volumetric heat generation (q’’’) to 
obtain a flux of 4 x 10° W/m*. Next the volumetric heat 
generation (q’’’) was constant which resulted in different heat 
fluxes from the chip. 

1. Constant Heat Flux (q’). 

The result of maintaining a constant heat flux is a 
chip that is a highly concentrated heat source when B=0.18 and 
a greatly distributed heat source when $B=0.81. In Figure 4.12 
the non-dimensional temperature profile is shown for B= 0.18 


to 0.81 at a value oti 
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Figure 4.12 Non-dimensional temperature for various size chips 
constant heat flux 

When the chip height results in a value of B=0.18, the 

chip is acting as a concentrated heat source with a volumetric 

heat generation q’” = 50 MW/m’ and the temperature of the chip 

is greatest. The larger the chip is to the rest of the 


package the lower the temperature of the overall package is. 


43 


A partial explanation for this is the amount of surface area 
exposed to convective cooling is smaller when B = 0.18 then 
when B = 0.81. * 

The thermal strain profile is shown if Figure 4.13 for 
the same value of ¢, = .564. The thermal strain profiles can 
be explained in much the same way the temperature profiles 


were. That is, when the chip is smallest the largest 


expansion of the package occurs. 
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Figure 4.13 Thermal strain profile for various chip heights 
constant heat flux 
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The temperature gradient 0, is shown in Figure 4.14 


for B values of 0.18 through 0.81. Figure 4.14 shows that 0, 
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Figure 4.14 6, for various chip heights at constant heat flux — 
increases with an increase in B, from ©, = 36 “C/mm at B = .18 
to 0, = 57 “C/mm at B = .81. When B = .18 the volumetric heat 
generation q’’ is largest at 50 MW/m’ and provides the least 
surface area exposed to convection. The thermal resistance 


for the substrate is attributable to the temperature gradient 
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developed in the substrate. At B = .18 the thermal resistance 


= = Ea ~ 
Reona (P=. 18) 16(.38)2 


is large, thereby not allowing much heat flux to conduct 


through the substrate, thereby keeping 0, small. 


At B = .81 the thermal resistance 
Rog (P=.81) =e 
rn - 16 (. 338) 


is very small, allowing for a much larger heat flux to pass 
through the substrate thereby increasing 0, Figure 4.15 


shows the temperature gradient 9, at varying 
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Figure 4.15 6, for various chip heights at constant heat flux 
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positions of ¢, for B of 0.18 through 0.81. These gradients 
are taken at the top of the substrate corresponding to a value 
et a= .47. At $, = .64 0, 1s approximately constant at 97 
°c/mm for all values of B. This can be attributed to the fact 
that no matter what thermal energy is generated by the chip it 
will tend to always flow through the solder and towards the 
upper right convective portion of the substrate. 

7 seene section of substrate between chips 
corresponding to ¢, = .75 and ¢, = .84, 9, increases slightly 
from 31 °C/mm to 52 °‘°C/mm and 14 °‘C/mm to 27 ‘°*C/mm 
respectively. This can be accounted for by looking at the 
substrate equivalent thermal resistances. When B = .81 the 
resistance is small in the Y direction as compared to the X 


direction as shown by, 


(sare i a Seen! a 


( Rees) y 


Therefore heat transfer would rather take place in the Y 
direction instead of the X direction. When B = .18, the 
thermal resistances are large in the Y direction and small in 


the X direction that is, 
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This change in thermal resistance now favors heat transfer in 
the X direction, which accounts for the imcrease 2n §@yyiaae 
smaller values of B. 
2. Constant Volumetric Heat Generation (q’’). 
This study shows the effect of using larger power 
chips in increasing sizes. The heat generation of each chip 
is proportional to chip volume, that is, heat generation 


increases with an increase in B. Figure 4.16 is the 
2 f0) 
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Figure 4.16 Non-dimensional temperature for various size chips 
constant volumetric power generation 
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normalized temperature profile for B = 0.18 to B = 0.81 ata 
position of ¢, = .564. 

The non dimensional temperature of the package follows 
the identical relationship as did volumetric heat generation 
to B. As B increases so does the non-dimensional temperature 
of the package. If the ratio of heat flux to non-dimensional 
temperature is examined, the increase in non-dimensional 
temperature is found to be non-linear. The following ratios 


q’/T, are given for various values of B, 





The increase in heat flux q’(W/m*) from B = .18 to B = .31 


is 12 KW/m* and from B = .6 to B = .81 the increase is 14 
KW/m*. Even though the increase in q’’ is 2 KW/m* more for the 
meer 68 the ratio of q’/T. remains constant. The factors 
that contribute to this are, an increase in the convection 
surface area of the chip and a decrease in the conduction 
thermal resistance of the substrate to allow for a greater 
heat transfer through the substrate. 

Figure 4.17 shows the thermal strain profile at a value of 
?, = .564. The chip does not have a large increase in thermal 
strain over the range of B’s while the solder has a very large 
increase in thermal strain. This goes back to Table 4.1, 


where the coefficient of thermal expansion for the Pb-Sn 
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solder is roughly eleven times greater than that of Silicon. 
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Figure 4.17 Thermal strain profile for various chip heights 
constant volumetric power 


Figure 4.18 shows the variation of the temperature 
gradient ©, at the top of the substrate at various values of 
py. At values of ¢, = 0.64 and ¢, = 0.43 the temperature 
gradient is largest and this corresponds to either side of the 
solder joint. At other values of ¢, the temperature gradient 


1s diminished greatly. At a value of B = 0.81 the temperature 
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gradients are large. They are almost twice as large as 
encountered in any other analysis. For this large a chip 
(large q’”) the convection cooling must be increased 
significantly to reduce this temperature gradient. For small 
chip sizes the temperature gradient is very small for any 


value of #¢,. 
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Figure 4.18 ®, for various chip heights at constant volumetric 
power 


spk 


Figure 4.19 shows the variation in 6, at the top of 
the substrate where the largest value of the temperature 
gradient within the depth of the substrate occurs. 

When 8B = 0.81 the chip has the greatest heat flux, 
q’'=72 KW/m* with a correspondingly large ®,. When B = 0.18 the 
chip has the lowest heat flux, q’= 16 KW/m* and 0, is all but 


non-existent. 
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Figure 4.19 8 for various chip heights at constant volumetric 
power 
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3. Overview Of Geometric Study 

For both the constant heat flux q” (W/m*) and the 
constant volumetric heat generation q’” (W/m*), the temperature 
gradient in the Y direction at the top of the substrate is 
greatest when the chip is largest (B = .81). The reasons as 
given previously are the interaction between the conductive 
thermal resistance, the convective thermal resistance, and the 
ratio of volumetric heat generation to available convection 
area. 

If the convective thermal resistance is small as 
compared to the conductive thermal resistance then the 
convective mode will dominate the heat transfer process as the 
heat transfer process is -governed by the path of least 


resistance. 


F. EFFECT OF CHIP ENCAPSULATION ON SYSTEM TEMPERATURE 

Some electronic packages are protected from moisture and 
the environment by a protective "encapsulating" coating. For 
this study the encapsulation material EME-1100-T.[Ref. 1] is 
0.1875 mm thick. The effect of encapsulation on the chip 
temperature and gradients developed within the package is 
studied. 

The following parameters will be held constant, volumetric 
heat generation (q’’’), convection cooling (h), and ambient 
temperature (T,). These values were the same as in section B 


for comparison purposes. The materials are a Silicon chip, a 
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Lead/Tin Solder, with all four substrates, and the 
encapsulation material. The thermal condttctivity Giga 
encapsulation material is k = 0.67 W/m°C. Within the other 
encapsulation materials the thermal conductivity (k) varied 
from 0.67 to 1.97 W/m C. 

Figure 4.20 is the non-dimensional temperature profile 
taken at .a value of ¢)— 900 55a. The entire package 
experiences less then a 5% temperature change, when low 
thermal conductivity substrates are used and less then 2% 
temperature change when high thermal conductivity substrates 
are used. The encapsulation material causes the surface 


temperature to decrease by a small amount. 
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Even though the encapsulation material has a low thermal 
conductivity its overall thermal resistance is low because its 


less then two tenths of a millimeter thick. 
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This 1s a very low thermal resistance accounting for the small 


change in temperature profile. Figure 4.21 shows the thermal 
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Figure 4.21 Encapsulated chip thermal strain profile 
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The thermal strain profile has the same negligible 
increases as the normalized temperature profile had except the 
encapsulation material has a very large coefficient of thermal 
expansion (a). 

In summary the effect of encapsulation on system 


temperature and temperature gradients are negligible. 
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V. STRESS ANALYSIS 


A. OVERVIEW OF STRESS ANALYSIS 

The stress analysis conducted is based on using different 
materials in the package and keeping all other factors 
constant. A silicon chip material was used for all cases, 
with two solders Pb-Sn and Au-Sn , and the four substrates 
Epoxy Fiberglass, Polyimide Fiberglass, Alumina and Aluminum 
Nitride. These eight sets of material combinations correspond 
to the thermal strain graphs, shown in Figures 4.3 and 4.4 in 
eaapter IV. 

Figure 4.3 corresponds to the material combinations of a 
Silicon chip, the Pb-Sn solder, and the four substrates listed 
above. Figure 4.4 corresponds to the material combinations of 
a silicon chip, the Au-Sn solder, and the four substrates 
listed above. 

The following other factors remain constant, q’’= 20 MW/m 
and h = 250 W/m°°C. All stresses given with the exception of 
the bending stresses at the centerline of the package are 
along the solder interfaces. The following stresses are 
calculated. First, the bending stresses (0,), along the 
package centerline. Three stresses are calculated along the 
solder-chip interface, the upper shear stress (Tt,), an upper 


bending stress (Op u), and an upper normal stress (o,,). Three 


218 | 


stresses are calculated along the solder-substrate interface, 
the lower shear stress (t,), a lower bending stress (Oo, .), and 
a lower normal stress (oO, ,)- 

Table 5.1 gives mechanical properties of electronic 
packaging materials, elastic modulus given in GPa. These 
properties are used in the calculation of stresses in 
accordance with the FEM development of chapter three. 


TABLE 5.1 


Elastic Modulus Poisson’s Ratio 


Silicon 13.03 
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B. STRESS ANALYSIS WITH (Pb-Sn) SOLDER 
1. Bending Stress Along Centerline, Oo, 

The bending stress at the centerline of the package is 
the lowest of all stresses presented and is shown in Figure 
5.1. The bending stress o,, through the chip vs of (Che s@aes 
of 1 MPa for all four substrates. Through the chip the 


bending stress is tensile for the low thermal conductivity 
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substrates (Epoxy and Polyimide Fiberglasses) and compressive 
for the high thermal conductivity substrates (Alumina and 


AlN). Through the substrates the bending stresses in the low 
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Figure 5.1 Centerline bending stresses o, 


thermal conductivity substrates 1s a maximum in compression at 
the top of the substrate (p, = -47) at -8 MPa and is a maximum 
in tension at the bottom of the substrate (p, = 0) at 6 MPa. 
The high thermal conductivity substrates have just the 


Opposite results, with a maximum in tension at the top of the 
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substrate “y= -47) at 2.5 Mpa and maximum in compression at 
the bottom of the substrate (py = 0) at -2.5 MPa. Thus the 
bending stresses in the low thermal conductivity substrates 
are about twice as large as the bending stresses in the high 
thermal conductivity substrates. The advantage with respect 
to bending stresses along the centerline of the chip is with 
the high thermal conductivity substrates (alumina and AIN). 
2. Upper Normal Stress, Fou 
The upper normal stress distribution along the chip- 
solder interface for all the substrates are compressive with 
a minimum occurring at 6 = .05 and maximum occurring at 6 = 
0.8 as shown in Figure 5.2. For the low thermal conductivity 
substrates the maximum upper normal stress for Epoxy 
Fiberglass is ~-49 MPa, and the maximum for Polyimide 
Fiberglass is -40 Mpa. Both curves show a decrease in the 
upper normal stress at the edges of the solder due to the fact 
that the edges at 6 = 0.0 and 6 = 1.0 are free surfaces. For 
the Pb-Sn group of materials the normal stress is the maximum 
stress examined for the low thermal conductivity substrates. 
The upper normal stresses for the high thermal 
conductivity substrates are about -16 MPa and are within 5% of 
each other, decreasing in magnitude slightly at the edges of 
the solder. This shows that the upper normal stresses in the 
high thermal conductivity substrates are about a third of the 


upper normal stress in the low thermal conductivity 


60 


substrates. We note that the normal stresses are fairly 
fiomtorm from 6 = «25 to b= .9. 

In contrast to the chip-solder interface, the normal 
stresses at the solder-substrate interface (i.e. the lower 
normal stresses) are very small on the order of KPa and in 
addition are tensile stresses for all substrate combinations. 
With respect to o,, the advantage is again seen with the high 


thermal conductivity substrates (alumina and Aln). 


10 
a CHIP 
a. 0 SOLDER 
= 
b° 
” -10 
VY) 
ul 
‘a 
y -20 
4 
$ 
oe -30 
oO 
2 
ac ~\_ _ _ POLYIMIDE FBGLS S 
“i -40 a ae ee 
a. 
~ EP/FBGLS 
-50 





NORMALIZED X LOCATION ALONG SOLDER 8s xb 
L 


Figure 5.2 Upper normal stresses a,, 
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3. Upper Bending Stress, o, , 

The upper bending stress o,, along the chip solder 
interface, is defined as the normal o, stress. However,-these 
stresses are presented along the interface, where they are 
greatest, as a function of 6 rather than across the Y 
direction through the solder. These bending stresses are 
compressive for all substrate combinations. For the low 
thermal conductivity substrates the maximum upper bending 
stresses occur at 6 = 0.05 and 0.95 with values of -16 MPa and 
-15 Mpa respectively. For the high thermal conductivity 
substrates maximum upper bending stresses occur at 6& = 0.05 
with -17 Mpa magnitude, and minimums of -8.5 MPa occur at 6 = 


0.8 as shown in Figure 5.3. Again we note that severe 
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Figure 5.3 Upper Bending Stress Cine 
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gradients occur close to the free edges, and with fairly 
uniform stresses along the center of the chip-solder 
interface. We note these bending stresses fall between the 
bending stresses at the centerline and normal stresses 
reported in the previous sections. With respect to these 
stresses, the advantage (i.e., lower stresses) is obtained for 
the low thermal conductivity substrates (epoxy and polyimide 
fiberglasses). 
4. Lower Bending Stress oa, |, 

The lower bending stresses are compressive in nature 
for all substrate combinations as shown in Figure 5.4. For 
the low thermal conductivity substrates, an immediate decrease 
from about -18 MPa to about -13 MPa, that is about a 30% 
change. Thereafter a gradual decrease to about -11 MPa, that 
is, about 15% change over the remaining domain. MThe lower 
bending stress for the high thermal conductivity substrates 
remains constant over the entire length of the solder joint at 
oo Pa. 

In comparison with the upper bending stresses the 
lower bending stresses are greater for the low thermal 
conductivity substrates. The advantage for the lower bending 
stresses goes with the high thermal conductivity substrates 


(alumina and AlN). 
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Figure 5.4 Lower bending stress 9, , 
5. Upper Shear Stress tr, 

The upper shear stresses at the chip/solder interface 
are positive at 6 = .05 and negative at 6 = .95 for all 
substrate combinations as shown in Figure 5.5. For the low 
thermal conductivity substrates the maximum shear stress is 
negative and occurs at 6 = .95 at -14 MPa and crosses zero 


stress at 6 = .15. For the high thermal conductivity 
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substrates the maximum shear stress is positive and occurs at 
6 = .05 at 11 MPa and crosses zero stress at 6 = .8. The 
upper shear stresses have large gradients at either edge of 
the chip-solder interface, but remain fairly constant through 
the majority of the domain. From 6 = .2 through 6 = .9 the 
upper shear stresses are less than 5 MPa in magnitude. 

These stresses fall in between the upper bending 
stresses and the bending stresses at the centerline. With 
respect to these stresses the advantage goes to the high 


thermal conductivity substrates (alumina and AI1N). 
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Figure 5.5 Upper shear stress tr, 
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6. Lower Shear Stress T, 

The lower shear stresses are opposite in sign from the 
upper shear stresses, in that the lower shear stresses at the 
solder/substrate interface are negative at 6 = .05 and are 
positive at 6 = .95 for all substrate combinations as shown 
in Figure 5.6. For the low thermal conductivity substrates 
the maximum shear stresses occur at 6 = .05 at about -16 MPa. 
This is comparable in magnitude with the upper shear stresses. 
For the high thermal conductivity substrates the maximum 
positive shear stress occurs at 6 = .95 at 16 MPa while the 
maximum negative shear stress occurs at 6 = .05 at -10 Mpa. 
This stress at 6 = .05 is twice the magnitude as compared to 
the upper shear stresses. Again we note that steep gradients 
occur at the edges of the solder-substrate interface and the 
low thermal conductivity substrates have a gradual decrease 
over the center of the domain. The t, stresses of Cicuaia 
thermal conductivity substrates cross zero at 6 = «9. 30a 
stresses for the high thermal conductivity substrates vary 
linearly with 6 between 6 = .2 and .9. These 1, Stteseee 
range from -5 MPa at 6 = .15 to 10 Mpa at & = .9 crossing zero 
stress at 6 = .45. We note these lower shear stresses are 
comparable in magnitude to the upper shear stresses. With 
respect to these stresses, there is no clear advantage for 


either set of substrates. 
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Figure 5.6 Lower shear stress T, 


C. STRESS ANALYSIS WITH (Au-Sn) SOLDER 
1. Bending Stress Along Centerline, oa, 

The bending stress at the centerline of the package 
remains the lowest of all stresses presented as shown in 
Figure 5.7. The bending stresses at the centerline (¢, = 0.0) 
for the packages with Au-Sn solder are identical to those with 


Pb-Sn solder. The centerline bending stress is independent of 
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which solder is used. With respect to these stresses, the 
advantage is obtained by the high thermal conductivity 


substrates (alumina and AlN). : 
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Figure 5.7 Centerline bending stress o, 
2. Upper Normal Stress, Per 

The upper normal stress distribution along the chip- 

solder interface for all the substrates are compressive in 


nature as shown in Figure 5.8. For the low thermal 


conductivity substrates the maximum for Epoxy Fiberglass is 
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-46 MPa and the maximum for Polyimide Fiberglass is -40 Mpa 


and they are constant over a range of .1 <6 < .9. Both upper 
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Figure 5.8 Upper normal stress 9g, , 
normal shear stress distributions decrease at the edges of the 
solder due to the fact that the stress at a free surface is 
zero. The normal stress is once again the maximum stress 
examined for the low thermal conductivity substrates. 

The high thermal conductivity substrates have 


increased by about 20% over the Pb-Sn material package. There 
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are sharp increases at 6 = .1 and a slight increase to a 
maximum of about -20 MPa at 6 = .9. 

Again we note that in contrast to the chip-solder 
interface, the normal stresses at the solder substrate 
interface are very small on the order of KPa and in addition 
are tensile stresses for all substrate combinations. With 
respect to these stresses, the advantage is obtained for the 
high thermal conductivity substrates (alumina and AI1N) and 
with respect to solders, the advantage is obtained for the Pb- 
Se solaen. 

3. Upper Bending Stress, 9, , 

The upper bending stress o,, along the chip-solder 
interface is defined as the normal o, stress as previously 
mentioned. The upper bending stresses are compressive for all 
substrate combinations as shown in Figure 5.9. For the low 
thermal conductivity substrates maximums occur at 6 = 0.05 and 
0.95 with values of -19 MPa and -22 Mpa respectively. For the 
high thermal conductivity substrates a maximum occurs at 6 = 
0.05 at -39 Mpa. This upper bending stress is twice as large 
for the high thermal conductivity substrates with Au-Sn solder 
then Pb-Sn solder. There is a very large gradient at 6 = 0.05 
where the upper bending stress decreases from -39 MPa to —-22 
MPa, that is, about a 45% decrease. For the high thermal 
conductivity substrates the upper bending stress continues to 


decrease at a gradual rate across the remaining domain, 6 = .2 
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to 6 = .9. With respect to these stresses, the advantage is 
obtained for the low thermal conductivity substrates (epoxy 
and polyimide fiberglasses) and with respect to solders, the 


advantage is obtained for the Pb-Sn solder. 
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Figure 5.9 Upper Bending Stress Oy u 
4. Lower Bending Stress, Opt 
The lower bending stresses are compressive for all 


substrate combinations as shown in Figure 5.10. For the low 


thermal conductivity substrates a maximum occurs at & = 0.05 


cal 


with a value of -39 MPa and a minimum at -7 MPa at 6& = 0.9. 
This is an increase of 200% over the same package with Pb-Sn 
solder. The lower bending stress for the high thermal 
conductivity substrates remain fairly constant over a range of 
.15 < 6 < .8 with a value of -38 MPa and is a maximum at 6 = 
.95 at -49 MPa. This is an increase of about 260% over the 
same package with Pb-Sn solder. The lower bending stresses 
are the highest stresses presented for the high thermal 
conductivity substrates. With respect to these stresses, the 
advantage is obtained for the low thermal conductivity 
substrates (epoxy and polyimide fiberglasses) and with respect 


to solders, the advantage is obtained for the Pb-Sn solder. 
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Figure 5.10 Lower bending stress a, , 
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5. Upper Shear Stress tT, 
The upper shear stresses at the chip/solder interface 
are positive at 6 = .05 and negative at 6 = .95 for all 


substrate combinations as shown in Figure 5.11. For the low 
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Figure 5.11 Upper shear stress tT, 

thermal conductivity substrates the maximum shear stress is 
negative and occurs at 6 = .95 with a value of -15 MPa and 
crosses zero stress at 6 = .15. For the high thermal 


conductivity substrates the maximum shear stress is positive 
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and occurs at 6 = .05 with a value of 16 MPa and crosses zero 
Stress at 6 = .8. The maximum upper shear stresses for both 
types of substrates show only a slight increase over the Pb-Sn 
solder packages. With respect to these stresses, there is no 
clear cut advantage for either set of substrates, nor either 
solder. 

6. Lower Shear Stress tT, 

The lower shear stresses are opposite in sign from the 
upper shear stresses, in that the lower shear stresses at the 
solder/substrate interface are negative at 6 = .05 and are 
positive at 6 = .95 for all substrate combinations as shown in 
Figure 5.12. For the low thermal conductivity substrates the 
maximum shear stresses occur at 6 = .05 with a value of -21 
MPa. This represents an increase of 130% over the same 
package with Pb-Sn solder. The lower shear stress then 
decreases to -10 Mpa at 6 = .1 and then slowly decreases to 
about zero at 6 = .95. For the high thermal conductivity 
substrates the maximum shear stresses occur at 6 = .95 with a 
value of 39 MPa. This represents an increase of 240% over the 
same package with Pb-Sn solder. At 6 = .05 the lower shear 
stress is negative at -32 MPa which represents an increase of 
110% over the same package with Pb-Sn solder. With respect to 
these stresses, the advantage is obtained for the low thermal 


conductivity substrates (epoxy and polyimide fiberglasses) and 
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with respect to solders, the advantage is obtained for the Pb- 


Sn solder. 
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Figure 5.12 Lower Shear Stress [, 


D. SUMMARY OF STRESS ANALYSIS 
The following table represents the maximum stresses (given 
in MPa) for each material combination studied. Various 


comparisons are possible. 


ifs 


TABLE 5.2 MAXIMUM STRESSES 
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For the Epoxy fiberglass substrate packages there is 
little difference in centerline bending o,, normal o,, and 
upper shear T, stresses and only a slightly larger upper 
bending o,,, and lower shear 1, stresses for Au-Sn versus Pb- 
Sn solder. For this substrate the lower bending o,, stresses 
are twice as great for Au-Sn versus Pb-Sn solder. The normal 
stresses are the largest stresses for this substrate for both 
Pb-Sn and Au-Sn solders. The most significant stresses with 
regard to deforming the solder joint are the lower bending and 
lower shear. The upper normal stress is not significant 
because it tends to compress a solder crack rather then 
propagate it. 

For the Polyimide fiberglass substrate packages there is 
little difference in centerline bending o,, normal o,, and 
upper shear 1T, stresses and only a slightly larger upper 
bending o,,, and lower shear 1, stresses for Au-Sn versus Pb- 
Sn solder. For this substrate the lower bending o,, stresses 
are twice as great for Au-Sn versus Pb-Sn solder. The normal 
stresses are the largest stresses for this substrate for both 
Pb-Sn and Au-Sn solders. The most significant stresses with 
regard to deforming the solder joint are the lower bending and 
lower shear. The upper normal stress is not significant 
because it tends to compress a solder crack rather then 
propogate it. The Polyimide fiberglass packages have stress 
distributions similiar to Epoxy fiberglass, but tend to be 


slightly lower in magnitude. 
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For the Alumina substrate only the centerline bending o, 
stresses are similiar for both solders. The upper normal o, 
and upper shear t, stresses are slightly larger for Aweenq 
versus Pb-Sn solders. The upper bending o,, stress is 225% 
larger, lower bending o,, stress is 300% larger, and the lower 
shear t, stress is 265% larger for the Au-Sn versus Pb-Sn 
solder. The upper bending stress is largest for the Pb-Sn 
solder while the lower bending stress is largest for the Au-Sn 
solder. The most significant stresses with regard to 
deforming the solder joint are the upper bending, lower 
bending and lower shear. 

For the AIN substrate only the centerline bending o, 
stresses are similiar for both solders. The upper normal o, 
and upper shear t, stresses are slightly larger for Au-Sn 
versus Pb-Sn solders. The upper bending Oo, stress is 225% 
larger, lower bending o,, stress is 305% larger, and the lower 
shear t, stress is 265% larger for the Au-Sn versus (Pea 
solder. The upper bending stress is largest for the Pb-Sn 
solder while the lower bending stress is largest for the Au-Sn 
solder. The most significant stresses with regard to 
deforming the solder joint are the upper bending, lower 
bending and lower shear. The AlN packages have stress 
distributions which are approximately less than 5% different 
from Alumina. 

Based on the stress results from ene investigation the 


following ranking of material combinations are presented based 
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on obtaining the lowest stresses in the package. The silicon 


chip is used in all packages. 
1. Pb-Sn with Alumina or Aln 
2. Pb-Sn with Polyimide fiberglass 
3. Pb-Sn with Epoxy fiberglass 
4. Au-Sn with Polyimide fiberglass 
5. Au-Sn with Epoxy fiberglass 
6. Au-Sn with Alumina or Aln 

Bending (o,) and Shear (1T) stresses at the solder- 
substrate interface can be significant. The lower inside 
corner of the solder-substrate interface (6 = 0.0) is the 
point where failure is most likely to occur due to bending and 
shear stresses. 

This investigator recommends the following topics for 
further research. Alter the thermal code for transient 
temperature response of the package. The result of this would 
be to track stress in time. In this way the changes in 
temperature and stress can be tracked in time after the 
electronic device has been turned on. Conduct a parametric 
study of a general tri-material package varying the materials 
coefficient of thermal expansion and modulus of elasticity to 


better understand the effect on stresses produced in the tri- 


material package. 


Ips, 


APPENDIX A. 


PROGRAM HEAT 
rrr ret er crrrcrcrrrrrrrcecreeeeeee ee eee eee eee eee eee eee eee eee eee eee ee ee ee | 


POISSON’S 


aalj 
ALPHA 
AREA 

bi 

bbij 
BIGA 
BIGF 
BDRYVAL 
BLOW 
cauchfi 
ci 
DEL1,2,3= 
DELMAX 
DIVSUB 
DIVX 
DIVY 
DIVXAIR 
DIVYAIR 
id 

GEN 
HLOW 
HUP 
HVERT 
IMAT 
ICOR 
KEL 
KINVS 
NBCLOW 
NBCUP 
NBCVERT 
NBDLEFT 


ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee 


NATREL 


THis: © see 
electronic packages. The governing differential equation is 


M. program solves for the Temperature distribution in 


equation. More specifically: 


When symmetry is envoked homogeneous Neumann boundary conditions 
are used at these locations. The other boundaries are Cauchy. 


Listing of variables: 


local nodal point values of BIGA matrix 


vector of element areas 

triangular coordinate property of each element 
contribution from Cauchy boundary to BIGA matrix 
global left hand matrix 

right hand side vector 

value of boundary point 

(BUP,BVERT) inverse of the absolute temp of the fluid 
contribution from Cauchy boundary to BIGF vector 
triangular coordinate property of each element 
difference between assumed and calc. surface temps 
maximum of dell or del2 or del3 

# of rows of substrate 

# of divisions in x direction for a grid density 

# of divisions in y direction for a grid density 

# of columns of air in grid 

# of rows of air(solder) in grid 

contribution of source term to BIGF vector 


if equal to 1 then internal heat generation is involved 


lower side heat transfer coefficient 
upper side heat transfer coefficient 
right side heat transfer coefficient 


+ © *£ & & FH HH HH HH He 


+ + + + HF HH HF HF HH HF HH HH KH HH 


identifies what material properties go with which element* 
establishes correspondence between local and global coor.* 


therm. conductivity of element 
kinematic viscosity of the fluid 
if = 1 ; then lower boundary is a Cauchy 


if = 1 ; then upper boundary is a Cauchy 
if = 1 ; then right boundary is a Cauchy 
if = 1 ; then left boundary is a Direchlet 
if = 1 ; then right boundary is a Direchlet 
if = 1 ; then upper boundary is a Direchlet 
if = 1 ; then lower boundary is a Direchlet 


number of columns for grid subroutine 
number of rows for grid subroutine 
# of substrate elements 
# of columns of air in grid 
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+e ee & & HF HE HE HF HF He FF 








# of air & solder elements 


RSOEL = 

Sar = volumetric heat generation of electronic device 

ALLOW = Raleigh # of the fluid at the lower surface 

ALUP = Raleigh # of the fluid at the upper surface 

ALVERT = Raleigh # of the fluid at the vertical surface 

40 = density of 4 different materials 

HOEL = density of element 

PHT = specific heat of 4 different materials 

PHTEL = specific heat of element 

AMB = ambient temperature 

=MP = solution vector of temperatures 

FLLOW = fluid temperature at lower surface 

FLUP = fluid temperature at upper surface 

FLVERT = fluid temperature at vertical surface 

SLOW = average surface temperature of lower surface 

5UP = average surface temperature of upper surface 

SVERT = average surface temperature of vertical surface 

JERCON = thermal conductivity of 4 different materials 
= x position of system nodal point 

208 = x position of new grid density 

= y position of system nodal point 

?0S = y position of new grid density 


hak keke KK KKK KKK KKK KK KKK KKK KKK KKK KK KK KKK KKK KKK KKK KKK KK KKK KKK KKK KKK 


kkk 
lek 
itk* 
itk&& 


|e 


itk & 


're 


| 


oe 
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INCLUDE ’HEATCOMN. FOR’ 


xx ensure the proper amount of workspace is allocated *** 
xx on problems with a large quantity of unknowns a ****xx 
xx FORTRAN STOP may be encountered. the amount of *****x 
** work space allocated in RWKSP and IWKIN must be **x*xx 
**x increased 

COMMON /WORKSP/RWKSP 

REAL RWKSP (900000) 

CALL IWKIN (900000) 


xx ensure the proper data files are being read ***kxxkkxr 
OPEN(7,FILE=’HEAT.OUT’ ,STATUS=’NEW’ ) 
OPEN(8,FILE=’TEST1.OUT’ ,STATUS='NEW’ ) 
OPEN(9,FILE=’TEST2.OUT’ ,STATUS='NEW’ ) 


ads given data and creates mesh 
CALL INPUT2 
CALL GRIDPLUS 


creates biga matrix and bigf vector 
CALL HEATMAT 


replaces Galerkin equation with direchlet bdry temperatures 
modifies equations so that heat flux at boundaries of material 
interfaces are equal 

CALL HEATMOD 


solves simultaneous system of equations 
CALL LSARG(NSNP,BIGA,NSNP,BIGF,1,TEMP) 


creates output file 
CALL OUTPUT 


END 
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* 
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* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
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* 
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SUBROUTINE INPUT2 


KK wm Kw KKK KRM KK KKK KKK KKK KKK KKK KKK KE KKK KKK KKK KEK KEK KKK KKEKKKKKKKAKKAKKKRKKK 
xke*k*k Reads in data to construct thermal grid 


INCLUDE ’HEATCOMN.FOR’ 
OPEN(11,FILE=’data8.DAT’ ,STATUS=’OLD’ ) 
OPEN(21,FILE=’TEST.OUT’ ,STATUS=’NEW’ ) 


READ(11,*) NCOL1,NCOL2,NROW - 
READ(11,*) (XPOS(I),I=1,NDIFX) 

READ(11,*) (YPOS(I),I=1,NDIFY) 

READ(11,*) (XINCOL(I),I=1,NDIFX-1) 

READ(11,*) (YINCRROW(I),I=1,NDIFY-1) 

READ(11,*) (THERCON(I),RHO(I),SPHT(I),I=1,NMATL) 
READ(11,*) DIVSUB,DIVXAIR, DIVYAIR, DIVXSOL 
READ(11,*) NBCLOW,NBCUP,NBCVERT, NBCSUB, TAMB 
READ(11,*) HLOW,HUP,HVERT,HSUB 

READ(11,*) GEN,QTPR 

READ(11,*) NBDUP,NBDRDEV,NBDUSUB,NBDRIT,NBDLOW 
READ(11,*) TUP,TRTDEV, TUPSUB, TRIGHT, TLOW 


END 
SUBROUTINE GRIDPLUS 


RH KKKHKKK KKK KKKKKKKKKKKKKAKKEKKEKKKKKKRKKKKKKRKKRKRKKKKA KKK KKK KKK KHK 


akkk overlays a grid on a square or rectangle with up to nine 
xxee different mesh densities in X and Y directions , sets up 
k*k*kk geometric and material properties for the Galerkin eq. 


INCLUDE ’HEATCOMN. FOR’ 


C open(21,file=’test.out’,status=’new’ ) 
DATA KEL, XVALUE, YVALUE /NEL*0,0,0/ 
RBMAX1 = NCOL1 + 1 
RBMAX2 = NCOL2 + 1 
RTMAX = NROW + 1 
XCOUNT = 1 
NN = ] 

RBMAX = RBMAX1 


KRHA KKKKKKKKKKKEKKKKKKKEKKKKKKAKKaAKAKKRK KKK KKK KK KKK KKK KKK KKK KKK KKAKRKK 
xxkk assigns X and Y coordinates to system Global nodal points 


DO 100 I = 1,NSNP 


X(I) = XVALUE 
Y(I) = YVALUE 


IF(XVALUE.LT.XPOS(2).AND.XVALUE.GE.XPOS(1))THEN 
XINCCOL = XINCOL(1) 

ELSEIF(XVALUE.LT.XPOS(3).AND.XVALUE.GE.XPOS(2)) THEN 
XINCCOL = XINCOL(2) 

ELSEIF(XVALUE.LT.XPOS(4).AND.XVALUE.GE.XPOS(3)) THEN 
XINCCOL = XINCOL(3) 

ELSEIF (XVALUE.LT.XPOS(5).AND.XVALUE.GE.XPOS( 4) ) THEN 
XINCCOL = XINCOL(4) 

ELSEIF(XVALUE.LT.XPOS(6).AND.XVALUE.GE.XPOS(5) )THEN 
XINCCOL = XINCOL(5) 
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ELSEIF(XVALUE.LT.XPOS(9).AND.XVALUE.GE.XPOS(8) ) THEN 
XINCCOL = XINCOL(8) 


ELSE 
XINCCOL = XINCOL(9) 
ENDIF 


XVALUE = XVALUE + XINCCOL 


IF( YVALUE.LT.YPOS(2).AND.YVALUE.GE.YPOS(1)) THEN 
YINCROW = YINCRROW(1) 

ELSEIF( YVALUE.LT.YPOS(3).AND.YVALUE.GE.YPOS(2) ) THEN 
YINCROW = YINCRROW( 2 

ELSEIF(YVALUE.LT. YPOS( 
YINCROW = YINCRROW( 

ELSEIF( YVALUE.LT. YPOS ( 
YINCROW = YINCRROW( 


4) .AND. YVALUE.GE.YPOS(3) )THEN 
3 
5 
4 
ELSEIF( YVALUE.LT.YPOS(6 
5 
7 
6 


»-AND.YVALUE.GE.YPOS(4)) THEN 


.AND.YVALUE.GE.YPOS(5)) THEN 
YINCROW = YINCRROW( 
ELSEIF(YVALUE.LT.YPOS( 
YINCROW = YINCRROW( 


.AND. YVALUE.GE.YPOS(6) ) THEN 


ELSEIF(YVALUE.LT.YPOS(8).AND.YVALUE.GE.YPOS(7)) THEN 
YINCROW = YINCRROW(7) 
ELSEIF( YVALUE.LT.YPOS(9).AND.YVALUE.GE.YPOS(8)) THEN 
YINCROW = YINCRROW( 8) 
ELSE 
YINCROW = YINCRROW(9) 
ENDIF 


begins numbering process of assembly after substrate 
IF(I.GT.((NCOL1+1)*(DIVSUB+1) ) ) THEN 
RBMAX = RBMAX2 
ENDIF 


IF( XCOUNT.GE.RBMAX) THEN 
XVALUE = 0.0 
YVALUE = YVALUE + YINCROW 
XCOUNT = 0.0 

ENDIF 


XCOUNT = XCOUNT + 1 
CONTINUE 


HORI KI HK KIKI KKK IK KI KK IK KKK IKI KIKI KK IKI KIKI KAKI KKK KKK KKK KK 
creates correspondence table for Local-Global nodal points 
DO 200 IROW = 1,DIVSUB + 1 
KK = (IROW-1)*(RBMAX1 ) 
Beeet0 ICOL # 1,NCOL1 


IEL = ((IROW-1)*(2*NCOL1) )+(2*ICOL)—-1 


web = IEL + 1 
| ICOR(IEL,1) = ICOL + KK 
| ICOR(IEL,2) = ICOL + NCOL1 + KK + 2 
ICOR(IEL,3) = ICOL + NCOL1 + KK + l 
ICOR(JEL,1) = ICOL + KK 
ICOR(JEL,2) = ICOL + KK + i 
| ICOR(JEL,3) = ICOL + NCOLL + KK + 2 
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210 CONTINUE 
200 CONTINUE 


*k*k* correspondence starting two rows above substrate 
KK = (RBMAX1*(DIVSUB+1) ) 
IEL = (NCOLI1*DIVSUB*2) + (NCOL2*2) + 1 
DO 220 IROW = 1,NROW-(DIVSUB+1) 
DO 230 ICGOL = 1PNECL2 


- 


JEL = IEL + 


ICOR(IEL,1) = ICOL + KK 

ICOR(IEL,2) = ICOL + NGOR2@yemee ee 

ICOR(IEL,3) = ICOL + NCOL2 + KK + l 

ICOR(JEL,1) = ICOL + KK 

ICOR(JEL,2) = ICOL + KK + l 

ICOR(JEL,3) = ICOL + NCOL2 + KK + 2 
Z 


LED = SPE Let 


230 CONTINUE 
KK = KK + RBMAX2 
220 CONTINUE 


KKK KKK KKK KK KKK KKK KKK KKK KKK KKK KK KKK KKK KKK KK KKK KK KAKA KKKKKKK KK KKK 
xkxkk assigns material properties to elements 


NSUBEL = NCOLI1 * DIVSUB * 2 
NAIREL = 2 * DIVXAIR 

NARSOEL = 2 * DIVYAIR * NCOL2 
LCOUNT = 1 

NCOUNT = l 

MCOUNT = 1 


DO 300 I = 1,NEL-1,2 
IF(I.LT.NSUBEL) THEN 
x*kk*x substrate (Single material) **** 
IMAT(I) 1 
TAGE (eee) 1 
x*kxk* substrate (multiple material) ***x 
ELSEIF(I.GT.NSUBEL.AND.I.LT. (NARSOEL+NSUBEL ) ) THEN 


IF (NCOUNT.LE.DIVXAIR) THEN 
kkk ajr kK 


IMAT (I) 2 
IMAT(I+1) 2 
ELSEIF(NCOUNT.LE.DIVXAIR+DIVXSOL ) THEN 
kk*k¥* solder **x** 


IMAT(I) = 3 
IMAT(I+1) = 3 
xkkk coating/cover *xxk 
ELSE 
IMAT(TI) = 5 
IMAT(I+1) = 5 
ENDIF 


NCOUNT = NCOUNT + 1 
IF (NCOUNT.GT.NCOL2 ) THEN 
NCOUNT = 1 
ENDIF 
ELSEIF(I.GT.NARSOEL+NSUBEL.AND. 
: I.LT.NEL~(NCOL2*2) ) THEN 
IF (MCOUNT.LE.NCOL2-1) THEN 
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device **xx 
IMAT(I) 
IMAT(I+1) 
coating/cover **** 
ELSE 
IMAT(I) 
IMAT(I+1) 
ENDIF 
MCOUNT = MCOUNT + 1 
IF(MCOUNT.GT.NCOL2 ) THEN 
MCOUNT = 1 


ll 
cr 


W 
63) 


5 


ENDIF 
ELSE 
coating/cover **** 
IMAT(TI) = 5 
IMAT(I+1) = 5 
ENDIF 
CONTINUE 
pO 310 I = 1,NEL 
II = IMAT(TI) 
KEL(I) = THERCON(II) 
RHOEL(I) = RHO(II) 


SPHTEL(I) = SPHT(II) 
IF(I.GT.(NSUBEL+NARSOEL) ) THEN 
QDOTEL(I) = QTPR 
ENDIF 
CONTINUE 


KKK KKK KKK KK KKK KKK KKK KKEKKKRKKKEKEKKEKKEKKKEKKEKKKEKRK KKK K KEK KK KK 


assigns boundary conditions to elements at the convective 
(CAUCHY) boundaries. 


DO 400 I = 2,2*NCOL1,2 


assigns Cauchy boundary to lower side of rectangle 
| IF(NBCLOW.EQ.1) THEN 
NBDRY(I) = l 
ENDIF 
CONTINUE 


DO 410 I = (NEL-(2*NCOL2))+1,NEL-1,2 


assigns Cauchy boundary to upper side of rectangle 
IF(NBCUP.EQ.1) THEN 
NBDRY(I) = 2 
ENDIF 

' CONTINUE 


DO 420 I = NSUBEL-( (NCOL1-NCOL2)*2)+1,NSUBEL-1,2 


assigns Cauchy boundary to upper tip of substrate 
| IF(NBCSUB.EQ.1) THEN 
NBDRY(I) = 4 
i ENDIF 
{ CONTINUE 
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DO 430 I = NSUBEL+(NCOL2*2),NEL,NCOL2*2 


k*x*k assigns Cauchy boundary to vertical side of rectangle 
IF(NBCVERT.EQ.1) THEN 
NBDRY(I) = 3 
ENDIF 
430 CONTINUE 


KH KKK KKH KK KKK KKK KK KK KK KKK KIKI KK KKK KIKI KKK KEKKE KK KKK KKK 


xxx assigns Direchlet boundary to system Global nodal points 
KKK 


kx*kk assigns Direchlet boundary to lower side of substrate 


DO 500 I = 1,RBMAX1 
IF(NBDLOW.EQ.1) THEN 


NDIRECH(I) = 1 
BDRYVAL(I) = TLOW 
ENDIF 


500 CONTINUE 


kxk*k assigns Direchlet boundary to right hand side of substrate 


DO 510 I = RBMAX1*2,RBMAX1*(DIVSUB+1),RBMAX1 
IF(NBDRIT.EQ.1) THEN 
NDIRECH(I) = 1 


BDRYVAL(I) = TRIGHT 
ENDIF 
510 CONTINUE 


x*kkk assigns Direchlet boundary to upper side of device 


DO 520 I = NSNP,NSNP-NCOL2,-1 
IF(NBDUP.EQ.1) THEN 
NDIRECH(1I) 
BDRYVAL(I) 
ENDIF 

520 CONTINUE 


1 
AI O))2 


x**k assigns Direchlet boundary to right side of device 
Ill = NSNP-(RBMAX2* (NROW-DIVSUB-1 ) ) 
DO 530 I = NSNP-RBMAX2, III,-RBMAX2 
IF(NBDRDEV.EQ.1)THEN 
NDIRECH(I) = 1 
BDRYVAL(I) = TRTIDEV 
ENDIF 
530 CONTINUE 


x*** assigns Direchlet boundary to upper side of substrate 


DO 540 I = (RBMAX1*DIVSUB)+RBMAX2, RBMAX1*(DIVSUB+1) 
IF(NBDUSUB.EQ.1) THEN 
NDIRECH(I) = 1 
BDRYVAL(I) = TUPSUB 
ENDIF 
540 CONTINUE 


END 
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SUBROUTINE HEATMAT 


te KK KK KH KKK KIKI KKK KKK KK KK KI KIRK KI KK KKK EK HKKKKKKKK KEKE KEKE 
calculates the BIGA matrix and BIGF vector that are the 
system of simultaneous equations to be solved; accounts 
for the main operator(Del”2), Cauchy boundary, and the 
forcing function to Poisson’s equation 


INCLUDE ’HEATCOMN.FOR’ 

DATA BIGA,BIGF / NDIM*0,NSNP*0 / 
mara bb21,bb22,bb23,bb33 /0.0,0.0,0.0,0.0/ 
DATA cauchfl,cauchf2,ff1 /0.0,0.0,0.0/ 


DO 10 IEL = 1,NEL 


Il = ICOR(IEL,1) 
I2 = ICOR(IEL,2) 
I3 = ICOR(IEL, 3) 
Bl = Y(I2) - Y(I3) 
B2 = Y(I3) - Y(I1) 
B3 = Y({I1) - Y(I2) 


= X(13) - X(1I2) 
@2-= X(1I1) - X(13) 
= X(I2) - X(I1) 


AREA(IEL) = (ABS(C1*B2-C2*Bl1))/2. 
write(8,*) iel,area(iel) 
aall = KEL(IEL)*(B1**2 + C1**2)/(-4.*AREA(IEL) ) 


aal2 = KEL(IEL)*(B1*B2 + C1*C2)/(-4.*AREA(IEL) ) 
aal3 = KEL(IEL)*(B1*B3 + ©C1*C3)/(-4.*AREA(IEL) ) 
aa21 = KEL(IEL)*(B2*Bl + C2*C1)/(-4.*AREA(IEL) ) 
aa22 = KEL(IEL)*(B2**2 + C2**2)/(-4.*AREA(IEL) ) 

aa23 = KEL(IEL)*(B2*B3 + C2*C3)/(-4.*AREA(IEL) ) 
aa3l = KEL(IEL)*(B3*Bl + C3*C1)/(-4.*AREA(IEL) ) 
aa32 = KEL(IEL)*(B3*B2 + ©C3*C2)/(-4.*AREA(IEL) ) 
aa33 = KEL(IEL)*(B3**2 + C3**2)/(-4.*AREA(IEL) ) 


KRkkkk kk aK KR KK KAKA KKRKKKA KKK RRR KKK KRKRKRKRKR KKK KKK KR KKK KK KKK KKK KR KKK 


calculation of Cauchy boundary to BIGA matrix and BIGF 
vector 


: IF(NBDRY(IEL) .NE.0) THEN 
boundary with local nodal points 1 - 2 on bottom side 
IF(NBDRY(IEL) .EQ.1) THEN 


|e) oh, = (HLOW*ABS(C3))/-3.0 

bb21 = (HLOW*ABS(C3))/-6.0 

cauchfl = (TAMB*HLOW*ABS(C3))/-2.0 
‘boundary with local nodal points 2 - 3 on upper side 


ELSEIF(NBDRY(IEL).EQ.2) THEN 


bb33 = (HUP*ABS(C1))/-3.0 
bb23 = (HUP*ABS(C1))/-6.0 
cauchf2 = (TAMB*HUP*ABS(C1))/-2.0 


‘boundary with local nodal points 2 - 3 on right side 
| ELSEIF(NBDRY(IEL).EQ.3) THEN 
bb33 = (HVERT*ABS(B1))/-3.0 
bb23 = (HVERT*ABS(B1))/-6.0 
a7 


(TAMB*HVERT*ABS(B1))/-2.0 


Gauchnez 


ELSE 
bb33 = (HSUB*ABS(Cl1))/-3.0 
bb23 = (HSUB*ABS(Cl1l))/-6.0 
cauchf2 = (TAMB*HSUB*ABS(C1))/-2.0 
ENDIF 
ENDIF 


KKK KH KKK KKK KKK KKK KK KKK KKK KKK KKK KK KKK KKK KKK KKK KKKKA KKK KKK KK 
kxx* the contribution of the Cauchy boundary 1s symmetric 


xkkkk bHbll=bb22 bbl2=bb21 bb22=bb33 bb23=bb32 
BIGA(I1,11) = BIGA(I1,I1) + aall bbZzZ 
BIGA(I1,1I2) = BIGA(I1,I1I2) + aal2 bb21 
BIGA(I1,13) = BIGA(I1,13) + aal3 
BIGA(I2,11) = BIGA(I2,I1) + aa2l bb21 
BIGA(I2,I2) = BIGA(I2,12) + aa22 bb22 + bb33 
BIGA(I2,13) = BIGA(I2,13) + aa23 bb23 
BIGA(I3,I1) = BIGA(I3,I11) + aa3l 
BIGA(I3,12) = BIGA(I3,I12) + aa32 bb23 
BIGA(I3,13) = BIGA(I3,13) + aa33 bb33 
bb22 = 0.0 
bb21 = 0.0 
Bbb235— 20 20 
bb33 = 0.0 


KKK KKRKK KK KKK KK KK KKK KKK KK KKK KKK KKK KKK KKKKKKKKKKAKKKAKKKKAKAKKK 
x*xxk used for Poisson’s equation (internal heat generation) 


IF(GEN.EQ.1.0) THEN 
IF( IMAT(IEL).EQ.4) THEN 
ffl = -AREA(IEL)*QTPREL(IEL)/3.0 
ENDIF 
ENDIF 


KKK KHKK KKK KKK KKK KK KKK KKK KK KKK KKK KKK KKK KKKKKK KAKA KK KKK KAKA KKKKK 
kxxk the two terms that comprise the BIGF matrix are from the 

x*kkk Cauchy boundary condition and the source term of the chip 
kxx*x a lumped approximation is used for both cases; cauchfl is 


xxxk for a 1-2 boundary edge cauchf2 is for a 2-3 boundary edge 
BIGF(I1) = BIGF(I1) + ff1 + cauchfl 
BIGF(I2) = BIGF(I2) + ff1 + cauchfl + cauchf2 
BIGF(I3) = BIGF(I3) + ff1 + cauchf2 
cauchfl = 0.0 
cauchtz #330. 0 
Be ry 0.0 

10 CONTINUE 


END 


SUBROUTINE HEATMOD 


KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KK KK KKK KKK KK KKKKKKKAKAKK KKK KKK KK KK 


x*xxx modifies BIGA matrix and BIGF vector if a system Global nodal 
kxxk point is a Direchlet boundary ie. replaces the Galerkin eq. 
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INCLUDE 'HEATCOMN.FOR 
modifies system of equations for Direchlet boundary points 


pO 10 I = 1,NSNP 


IF (NDIRECH(I) .EQ. 1)THEN 
BIGA(I,I) = 1.0 


po 20 J = 1,I-1 
BIGA(I,J) = 0.0 
CONTINUE 


DO 30 JJ = I+1,NSNP 
BIGA(I,JJ) = 0.0 


CONTINUE 
BIGF(I) = BDRYVAL(I) 
ENDIF 
CONTINUE 


END 
SUBROUTINE OUTPUT 


RKEKEKKKKKKEKK KKK KKK K KKK HK KKK KKK KK KKEK KKK KHER KKK KKK KKEKKAKKKKKKKKAK KKK 


creates output file of temperature distribution and a data file 
for graphing the temperature contours 


INCLUDE ’HEATCOMN. FOR’ 
OPEN(10, ee see DAT’ ,STATUS=’ NEW’ ) 
write(7,*) Remeeeaicares through solder’ 
write(7,*) Per tenet 8). eo, eemp 235) ,' 36) ,temp(38),'53’, temp(53 ) 
write(7,*) ‘68’, temp(68),’83’, temp(83),’93', temp(93), 
: ’103’,temp(103) 
write(7,*) '113’,temp(113),’'123’,temp(123),'133’,temp(133) 


energy balance 


gqconv = 0.0 
do 200 i= 126,134 
qconv=qconv+(hup*(((temp(I)+temp(1I+1))/2)-tamb) *(X(I+1)-x(I))) 


continue 

do 210 i= 1,14 

qconv=qconv+(hlow*(((temp(I)+temp(I+1))/2)- —-tamb) *(X(I+1)-x(I))) 
continue 

Zeee20 i = 70,74 

qconv= qconv+(hsub*(((temp(I)+temp(1+1))/2)- -tamb) *(X(I+1)-x(I))) 
continue 

Boe230 i = §85,125,10 

Aconvegconv+(hvert*(( ( temp(1)+temp(I+10) ) /2)- tamb ) * 


(y(I+10)-y(I))) 
‘continue 
qconv=qconv+(hvert*(((temp(70)+temp(85))/2)-tamb) *(y(85)-y(70))) 


WRITE(7, *) 
write(7,*) ’ Heat loss to convection = ',qconv,’W/M’ 


POWER = QPRT*0.0025*0.002 
Bec l as ,»*) ‘The power output of the device =’,power,’W/M’ 
ND 
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APPENDIX B. 


PROGRAM FILTER 
REAL X(135),Y(135),T(135) ,TEMPT(200)) [eMPee og 


OPEN(12,FILE=’bh3.out’ ,STATUS=’OLD’ ) 
OPEN(13,FILE=’bh3f£.DAT’ ,STATUS=’NEW’ ) 


READ(12,*) (X(1), Y( 1) 30 Gee 


xxx*k* calculates intermediate temperatures ***xxxx 


TA = (T(61)-T(46))*.75 + 1T(46) 

TB = (T(61)+T(46))/2.0 

TC @? (T(62)-T(47))*. 75 F Tea 

TD = (T(66)+T(47))/2.0 

TE = (T(63)-T(48))*.75 + 1T( 48) 

TF = (T(63)+T(48))/2.0 

TG = (7T(64)-T(49))*.75 + T(49) 

TH = (1T(64)+T(49))/2.0 

TI = (T(65)-T(64))*.75 + 1T(64) 

TL = (T(50)-T(49))*.75 + 1T(49) 

TJ = (TI-TL)*.75 + 2 

TK gee (ee Db 72 

TM = (T(65)-T(50))*.75 + T(50) 

TN = (7T(65)+T(50))/2.0 

TO ty (T(65) 466 yen Oo 

TR = (T(51)-T(50))*.75 + 1T(50) 

TP = (TO-TR)*.75 + TR 

TQ = (TO+TR)/2.0 

TS = (1T(65)+T(66))/2. 

TV = (T(50)4+T(51))/2. 

TT = (TS=TV)*.75°4 =i 

TU = (TS+TV)/2.0 

TW = (7T(66)-T(51))*.75 + T(51) 

TX = (7T(66)4+T(51))/2.0 

TY = (T(67)-T(52))*.75 + T( 52) 

TZ = (T(67)4+T(52))/2.0 

TAA = (T(68)-T(53))*.75 + T(53) 
TAB = (T(68)+T(53))/2.0 

TAC = (1T(69)-T(54))*.75 + 17(54) 
TAD = (1T(69)+T(54))/2.0 

TAE = (T(69)+T(70))/2.0 

TAH = (7T(54)4+T(55))/2.0 

TAF = (TAE-TAH)*.75 + TAH 

TAG = (TAE+TAH) /2.0 

TAT tare( T(69)—T(70) )*. 75 eee 0) 
TAL = (1T(54)-7T(55))*.75 + 1(54) 
TAJ *# ( TAI=TAD)% 75a a 

TAK = (TAI+TAL) /2.0 

TAM = (T(70)-T(55))*.75 + T(55) 
TAN = (1T(70)+T(55))/2.0 

TAO ge o(T(70)-T( 71) )* 5) era 
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TAR = (T(55)-7T(56))*.75 + 1T(55) 
TAP = (TAO-TAR)*.75 + TAR 

TAQ = (TAO+TAR)/2.0 

TAS = (T(71)-T(56))*.75 + 1T(56) 
TAT = (T(71)+T(56))/2.0 

TAU = (1T(75)-T(60))*.75 + 1T(60) 
TAV = (1T(75)+T(60))/2.0 

TAW = (T(35)-T(34))*.75 + T(34) 
TAX = (T(36)-T(35))*.75 + 1T(35) 
TAY = (T(36)+T(35))/2.0 

TAZ = (T(39)+T(40))/2.0 

TBA = (T(39)-T(40))*.75 + T(40) 
TBB = (1T(40)-T(41))*.75 + 1T(41) 
TBC = (T(5)-T(4))*.75 + T(4) 
TBD = (T(6)-T(5))*.75 + T(5) 
TBE = (T(6)+T(5))/2.0 

TBF = (T(9)+T(10))/2.0 

TBG = (T(9)-T(10))*.75 + T(10) 
TBH = (T(10)-T(11))*.75 + T(11) 


calculates average temperatures for rectangular element FEM program 
yeld.for (chip and solder) 


TEMPS1 = (7T(95)+7T(85)+T(96)+T(86))/4.0 
DO 100 I = 1,115 

TEMPT(I) = TEMPS1 

TEMPB(I) = TEMPS1 
CONTINUE 


calculates average temperatures for rectangular element FEM program 
yveld.for (substrate) 


TEMPT(116) 
/TEMPB(116) 


(T(61)4+T(61))/2.0 
(TA+TC)/2.0 


TEMPT(117) 
-TEMPB(117) 


(T(62)+T(63))/2.0 
(TC+TD)/2.0 


‘TEMPT(118) 
TEMPB(118) 


= (T(63)+T(64))/2.0 
= (TE+TG)/2.0 


TEMPT(119) 
TEMPB(119) 


‘TEMPT(120) 
TEMPB(120) 


TEMPT(121) 
TEMPB(121) 


TEMPT( 122) 
TEMPB(122) 


TEMPT(123) 
TEMPB(123) 


-TEMPT(124) 
TEMPB(124) 


(T(64)4+TI)/2.0 
(TG+TJ)/2.0 


= (T(65)4+TI)/2.0 
= (TJ+TM)/2.0 


il 


(T(65)4+T0)/2.0 
(TM+TP)/2.0 


(TO+TS)/2.0 
(TP+TT) /2.0 


(T(66)+TS)/2.0 
(TT+TW)/2.0 


(T(66)4+T(67))/2.0 
(TW+TY)/2.0 
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TEMPT(125) 
TEMPB(125) 


TEMPT(126) 
TEMPB(126) 


TEMPE TZ 7) 
TEMPB(127) 


TEMPT(128) 
TEMPB(128) 


TEMPT(129) 
TEMPB(129) 


TEMPT(130) 
TEMPB(130) 


TEMPT(131) 
TEMPB(131) 


TEMPT( 132) 
TEMPB(132) 


TEMPT(133) 
TEMPB(133) 


TEMPT(134) 
TEMPB(134) 


TEMPT(135) 
TEMPB(135) 


TEMPT(136) 
TEMPB( 136) 


TEMPT( 137) 
TEMPB(137)} 


TEMPT( 138) 
TEMPB(138) 


TEMPT(139) 
TEMPB(139) 


TEMPT(140) 
TEMPB(140) 


TEMPT(141) 
TEMPB(141) 


TEMPT(142) 
TEMPB(142) 


TEMPT(143) 
TEMPB(143) 


TEMPT(144) 
TEMPB(144) 


(T(67)+TC6S 20 
(TY+TAA) /2.0 


(T(68)+T(69))/2.0 
(TAA+TAC) /2.0 


(T(69)+TAE) /2.0 
(TAC+TAF) /2.0 


(TAE+TAI)/2.0 
(TAF+TAJ)/2.0 


(T(70)+TAL Varo 
(TAJ+TAM) /2.0 


(T(70)+TAO) /2.0 
(TAM+TAP) /2.0 


(T(71)+TAO)/2.0 
(TAP+TAS) /2.0 


(Tig dict aSaay 2 20 
(TAS+TAU)/2.0 


(TA+TC)/2. 


= ( TE+ TE) /2.. 
= (TD+TF)/2. 


ot 


= (TAA+TAC) /2. 
= (TAB+TAD)/2. 


0 
(TB+TD) /2.0 
0 
0 


(TE+TG) /2.0 
(TF+TH)/2.0 


(TG+TJ)/2. 
(TH+TK) /2. 


(TI+TM) /2.0 
(TK+TN) /2.0 


(TM+TP)/2.0 
(TN+TQ) /2.0 


(TP+TT)/2.0 
(TQ+TU) /2.0 
0 
0 


(TT+TW) /2. 
(TU+TX) /2. 


(TW+TY)/2.0 
(TX+TZ)/2.0 


(TY+TAA) /2.0 
(TZ+TAB) /2.0 


(TAC+TAF) /2. 
(TADENG) /2- 


> a) oo >) 


oz 


TEMPT(145) 
TEMPB(145) 


TEMPT(146) 
TEMPB(146) 


TEMPT(147) 
TEMPB(147) 


TEMPT(148) 
TEMPB(148) 


TEMPT(149) 
TEMPB(149) 


TEMPT(150) 
TEMPB(150) 


TEMPT(151) 
TEMPB(151) 


TEMPT(152) 
TEMPB(152) 


TEMPT(153) 
TEMPB(153) 


TEMPT(154) 
TEMPB(154) 


TEMPT(155) 
TEMPB(155) 


TEMPT(156) 
TEMPB(156) 


TEMPT(157) 
TEMPB(157) 


TEMPT(158) 
TEMPB(158) 


TEMPT(159) 
TEMPB(159) 


TEMPT(160) 
TEMPB(160) 


TEMPT(161) 
TEMPB(161) 


TEMPT(162) 
TEMPB(162) 


TEMPT(163) 
TEMPB(163) 


TEMPT(164) 
TEMPB(164) 


= (TAS+TAU)/2. 
= (TAT+TAV)/2. 


(TAF+TAW) /2. 
(TAG+TAK)/2. 


(TAJ+TAM)/2. 
(TAK+TAN) /2. 


Oo © a> i a) 


(TAM+TAP) /2.0 
(TAN+TAQ)/2.0 


(TAP+TAS)/2. 
(TAQ+TAT) /2. 


Oo © Oo © 


(TB+TD)/2.0 
(T(46)4+T(47))/2.0 


(TD+TF)/2.0 
(T(47)4+T(48))/2.0 


(TF+TH)/2.0 
(T(48)+T(49))/2.0 


(TH4+TK)/2.0 
(TL+T(49))/2.0 


(TK+TL)/2.0 
(TN+T(50))/2.0 


= (TN+TQ)/2.0 
= (TR+T(50))/2.0 


(TQ+TU) /2.0 
(TR+TV) /2.0 


(TU+TX)/2.0 
EVs T Goan) /Z.0 


(TX+TZ)/2.0 
eo eres 2.) 7 2 « 0 


(TZ+TAB)/2.0 
(T(52)+T(53))/2.0 


(TAB+TAD)/2.0 
(T(53)+T(54))/2.0 


= (TAD+TAG) /2.0 
= (T(54)+TAH)/2.0 


(TAG+TAK)/2.0 
(TAH+TAL) /2.0 


(TAK+TAN)/2.0 
(TAL+T(55))/2.0 


(TAN+TAQ) /2.0 
(TAR+T(55))/2.0 
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TEMPT(165) 
TEMPB(165) 


TEMPT( 166) 
TEMPB(166) 


TEMPT( 167) 
TEMPB(167) 


TEMPT(168) 
TEMPB(168) 


TEMPT(169) 
TEMPB(169) 


TEMPT(170) 
TEMPB(170) 


TEMPT(171) 
TEMeB( 171) 


TEMPT(172) 
TEMPE? 2) 


TEMPT(173) 
TEMPB(173) 


TEMPT(174) 
TEMPB(174) 


TEMPT(175) 
TEMPB(175) 


TEMPE TD (AG) 
TEMPB(176) 


TEMS Tele) 
TEMPB(177) 


TEMPT(178) 
TEMPB(178) 


TEMPT(179) 
TEMPB(179) 


TEMPT(180) 
TEMPB(180) 


TEMPT(181) 
TEMPB(181 ) 


TEMPT(182) 
TEMPB(182) 


TEMPT( 183) 
TEMPB(183) 


TEMPT(184) 
TEMPB(184) 


(TAQ+TAT) 72.0 
(TAR+T(56))/2.0 


(TAT+TAV)/2.0 


(T(56)+T(60) )/2. 


(T(49)4TL)/2. 
(T(34)+TAW) /2.0 


4 

3 
50)+TL)/2.0 
35)+TAW)/2.0 


a 
ay 


(T(50)4+TR)/2.0 
(T(35)+TAX) /2.0 


(TR+TV) /2.0 
(TAX+TAY)/2.0 


)+TV)/2.0 


(T-5 
(T(36)+TAY)/2.0 


2 
2 
(TAH+TAL) /2.0 
(TAZ+TBA) /2.0 


(T(55)+TAL) /2.0 
(T(40)+TBA)/2.0 


(T(55)4+TAR) /2.0 
(T(40)+TBB)/2.0 
0 
0 


(T(56)+TAR)/2. 
(T(41)+TBB7 2. 


on ae) Oo © 


oOo & Oo © 


Oo © 


Oo © 
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TEMPT(185) 
TEMPB(185) 


TEMPT( 186) 
TEMPB( 186) 


TEMPT(187) 
TEMPB(187) 


TEMPT(188) 
TEMPB(188) 


TEMPT(189) 
TEMPB(189) 


TEMPT(190) 
TEMPB(190) 


TEMPT(191) 
TEMPB(191) 


TEMPT(192) 
TEMPB(192) 


TEMPT(193) 
TEMPB(193) 


TEMPT(194) 
TEMPB(194) 


TEMPT(195) 
TEMPB(195) 


TEMPT(196) 
TEMPB(196) 


TEMPT(197) 
TEMPB(197) 


TEMPT(198) 
: TEMPB(198 ) 


TEMPT(199) 
TEMPB(199) 


(iis a5 (33))72.0 


(T(2)+7T(3))/2.0 


(T(33)+T(34))/2. 


(T(3)+T(4))/2.0 


(T(34)+TAW) /2.0 
(T(4)+TBC)/2.0 


(T(35)+TAW) /2.0 
(T(5)+TBC)/2.0 


(T(35)+TAX)/2.0 
(T(5)4+TBD)/2.0 


X+TAY)/2.0 
D+TBE) /2.0 


w YP 


(T(36)+TAY)/2.0 
(T(6)+TBE)/2.0 


36)4+T(37))/2. 


(T( 
(T(6)+T(7))/2.0 


CPs abies )) 722 


(T(7)+T(8))/2.0 


aH 


( 
(T(8)4+T(9))/2.0 
(T(39)+TAZ)/2.0 
(T(9)+TBF)/2.0 


(TAZ+TBA) /2.0 
(TBF+TBG) /2.0 


40)+TBA)/2.0 
10)+TBG)/2.0 


( 

(T(10)+TBH)/2.0 

41)+TBB)/2.0 
1 


( T( 
(T(11)+TBH)/2.0 


(oie ce Wes 7 2. 
( 


0 

0 
T(40)+TBB)/2.0. 

(10 


TEMPT( 200) 
TEMPB(200) 


"N 
| 


41)+T(45))/2.0 
(T(11)4+T(15))/2.0 


‘use 1f£ printout is desired 
DO 900 I = 1,200 

WRITE(13,*) TEMPT(I),’ ‘°,TEMPB(I) 

) CONTINUE 


END 
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